Mapping potato late blight from UAV-based multispectral imagery
Part 1
Jorge Luis Rodríguez
email: jorodriguezga@unal.edu.co
Introduction¶
This work proposes a method to process high resolution multispectral images for detection and mapping late blight in potato crops. It is based on three main aspects: (i) radiometric and geometric correction of the raw images; (ii) morphological operations for background removing; (iii) model trainig and supervised classification stage using a random forest classifier; and (iv) implementation of trained model in a new dataset.
Data and methods¶
The study area is a 1920 sq.m. potato field located in Subachoque, Colombia (see figure). This field is part of an experimental plot designed to evaluate the potato Diacol capiro variety response to different nutrient treatments. As the weather conditions favoured the appearance of late blight disease in several plots, the project had the chance to monitor the crop disease development. The experimental plot was inspected every seven days during the entire crop life span (i.e. 120 days from planting to maturity).
Mosaic generated extend over 3.2 Ha, which include the experimental plot and a bigger area where a different variety of potato was located, being this variety a yellow potato, this means, the tuber root colour it produces is yellow. In contrast, the potato from the experimental plot was from one variety of white potato. The experimental plot area was clipped by using a reference polygon (red polygon).
Study area:¶
High resolution multispectral images were acquired at 40 m altitude above the ground surface at 11:00 am local time (GMT-5). Each multispectral image acquired by the MicaSense camera had five bands as described in the next table. In this work camera bands 4 and 5 have been reset according to the following order: Blue (B); Green (G); Red (R); Red edge (RE) and Near infrared (NIR).
| # Banda | Nombre | Centro de Banda (nm) | Ancho Banda FWHM (nm) |
|---|---|---|---|
| 1 | Azul | 475 | 20 |
| 2 | Verde | 560 | 20 |
| 3 | Rojo | 668 | 10 |
| 4 | Infrarrojo cercano |
840 | 40 |
| 5 | Borde Rojo | 717 | 10 |
Experimental crop area consisted of an array of 18 blocks of 12 m $\times$ 8 m inside of a field of 77 m $\times$ 24 m. Each block had 9 rows arranged along the long side of the field with an area of 216 sq.m., the space between rows was 1m. Each row contained approximately 30 seed tubers of the same variety.
(a) May 12, 2018
(b) May 26, 2018
The layout of the experimental crop. The yellow lines indicate the division of the field into 18 experimental blocks. (a) May 12, 2018; (b) May 26, 2018.
In [1]:
%matplotlib inline
import time
import os, shutil
from matplotlib import pyplot as plt
from matplotlib import colors
from IPython.display import Image
import numpy as np
import cv2
import math
import time
from osgeo import gdal
from skimage import img_as_ubyte
import scipy
import pandas as pd
from skimage.filters import threshold_otsu
import skimage.io as io
from sklearn.ensemble import AdaBoostClassifier, RandomForestClassifier, GradientBoostingClassifier, ExtraTreesClassifier
from sklearn.externals import joblib
from sklearn.metrics import plot_roc_curve, classification_report, multilabel_confusion_matrix
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC, LinearSVC
from sklearn.neighbors import KNeighborsClassifier
In [2]:
ds = gdal.Open('Raster/Correccion_reflectancia/MR_CE_20180512_Subachoque.tif')
In [3]:
# asignacion cada banda
AZUL = ds.GetRasterBand(1).ReadAsArray()
VERDE = ds.GetRasterBand(2).ReadAsArray()
ROJO = ds.GetRasterBand(3).ReadAsArray()
REDEDGE = ds.GetRasterBand(4).ReadAsArray()
NIR = ds.GetRasterBand(5).ReadAsArray()
srs = ds.GetProjectionRef()
geo_transform = ds.GetGeoTransform()
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(321) ,plt.imshow(AZUL, cmap='gray'),plt.title('Blue band')
plt.subplot(322) ,plt.imshow(VERDE, cmap='gray'),plt.title('Green band')
plt.subplot(323) ,plt.imshow(ROJO, cmap='gray'),plt.title('Red band')
plt.subplot(324) ,plt.imshow(REDEDGE, cmap='gray'),plt.title('Red edge band')
plt.subplot(325) ,plt.imshow(NIR, cmap='gray'),plt.title('Infrared band')
plt.show()
Thresholding¶
A thresholding technique reduces a grey-level image into an image where objects and background are represented with two levels: a binary image (Glasbey, 1993). Because of the difference in reflectance between the soil and the potato plants, both near infrared and red edge bands allowed separation of the potato plants from the soil.
The first step was to analyse multispectral bands of the orthomosaics by plotting their histograms. The next figure shows the histogram for each band in the multispectral orthomosaic for the two categories under study as well as ground reference images of the two categories. Green, Rededge and Near infrared bands have right-skewed bimodal histograms in the first category. Bands with high weed presence had a right-skewed histograms without a clear valley. Thresholding method used to separate vegetation was OTSU's thresholding method (Otsu, 1979).
In [4]:
hist_b1 = np.histogram(AZUL[AZUL!=0].ravel(),100,[0,1])[0]
hist_b2 = np.histogram(VERDE[VERDE!=0].ravel(),100,[0,1])[0]
hist_b3 = np.histogram(ROJO[ROJO!=0].ravel(),100,[0,1])[0]
hist_b4 = np.histogram(REDEDGE[REDEDGE!=0].ravel(),100,[0,1])[0]
hist_b5 = np.histogram(NIR[NIR!=0].ravel(),100,[0,1])[0]
hist_range = np.arange(0,1,0.01)
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231), plt.bar(hist_range,hist_b1,0.01), plt.title("Blue")
plt.subplot(232), plt.bar(hist_range,hist_b2,0.01),plt.title('Green')
plt.subplot(233), plt.bar(hist_range,hist_b3,0.01),plt.title('Red')
plt.subplot(234), plt.bar(hist_range,hist_b4,0.01),plt.title('Rededge')
plt.subplot(235), plt.bar(hist_range,hist_b5,0.01),plt.title('Near infrared')
Out[4]:
In [5]:
image=NIR
hist_image = image.ravel()
hist = np.histogram(image[image!=0].ravel(),100,[0,1])[0]
thresh = threshold_otsu(image)
binary = image > thresh
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(131) ,plt.imshow(NIR, cmap='gray'),plt.title('Near infrared')
plt.subplot(132) ,plt.imshow(binary,cmap='gray'),plt.title('Binary image')
plt.show()
plt.show()
Masking image¶
To delete background from original multispectral image we multiply each band for the binary image resulting from thresholding step using the next equation:
\begin{equation} \label{Eq3} C_{i}(x,y) = A_{i}(x,y)B(x,y), \end{equation}where $C_{i}(x,y)$ is the resulting pixel at position $(x,y)$ for $i$ band without background, $A_{i}(x,y)$ is the original pixel at position $(x,y)$ for $i$ band of the multispectral image and $B(x,y)$ is the binary image pixel at position $(x,y)$ created in the thresholding process.
In [6]:
mask = np.array(binary, dtype='uint8')
filtered_mask = cv2.medianBlur(mask,3)
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(131) ,plt.imshow(NIR, cmap='gray'),plt.title('Banda NIR')
plt.subplot(132) ,plt.imshow(mask,cmap='gray'),plt.title('Mascara')
plt.subplot(133) ,plt.imshow(filtered_mask,cmap='gray'),plt.title('Mascara con filtro de mediana')
plt.show()
plt.show()
Here we can see how looks any band from the multispectral image after masking process. It can be seen how the ground surface does not appear any more.
In [7]:
azul_seg = AZUL*filtered_mask
verde_seg = VERDE*filtered_mask
rojo_seg = ROJO*filtered_mask
rededge_seg = REDEDGE*filtered_mask
nir_seg = NIR*filtered_mask
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(321) ,plt.imshow(azul_seg, cmap='gray'),plt.title('Blue')
plt.subplot(322) ,plt.imshow(verde_seg, cmap='gray'),plt.title('Green')
plt.subplot(323) ,plt.imshow(rojo_seg, cmap='gray'),plt.title('Red')
plt.subplot(324) ,plt.imshow(rededge_seg, cmap='gray'),plt.title('Red edge')
plt.subplot(325) ,plt.imshow(nir_seg, cmap='gray'),plt.title('Near infrared')
plt.show()
We can now create a false colour composition to improve the visualization of the background removing step. This was performed for visualisation purposes only:
In [8]:
# Contrast stretch function for enhance the bands
# Only for visualization purposes
def enhance_band(band, min_in= 0.1, max_in= 0.2):
min_setup = 0.0
max_setup = 1.0
enhanced = (band-min_in)*(((max_setup-min_setup)/(max_in-min_in))+min_setup)
return enhanced
In [9]:
NIR_IN = enhance_band(NIR, min_in= 0.0, max_in= (NIR.max()-0.12))
ROJO_IN = enhance_band(ROJO, min_in= 0.0, max_in= (ROJO.max()-0.06))
VERDE_IN = enhance_band(VERDE, min_in= 0.0, max_in= VERDE.max())
rgb = np.stack([NIR_IN,ROJO_IN,VERDE_IN], axis=2) #axis =2 so its shape is (M,N,3) otherwise doesn't work with plt. (M, N, 3): an image with RGB values (0-1 float or 0-255 int)
mask = cv2.cvtColor(filtered_mask, cv2.COLOR_GRAY2BGR) # So we can apply the mask
plantsImage = rgb*mask
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(121) ,plt.imshow(rgb),plt.title("Original image")
plt.subplot(122) ,plt.imshow(plantsImage),plt.title("Masked image")
plt.show()
plt.figure(2, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(121) ,plt.imshow(rgb[1450:1700,1100:1300]),plt.title("Original image")
plt.subplot(122) ,plt.imshow(plantsImage[1450:1700,1100:1300]),plt.title("Masked image")
plt.show()
Segmented bands to files¶
Now we proceed to save each of the segmented bands:
In [10]:
bands_out = np.stack([azul_seg,
verde_seg,
rojo_seg,
rededge_seg,
nir_seg],
axis=2)
band_names = {'M_2018-05-12_Corte_Seg_B1.tif':0,
'M_2018-05-12_Corte_Seg_B2.tif':1,
'M_2018-05-12_Corte_Seg_B3.tif':2,
'M_2018-05-12_Corte_Seg_B4.tif':3,
'M_2018-05-12_Corte_Seg_B5.tif':4
}
dir_out = 'Raster/Segmented/'
for band_name in band_names:
filename_output = dir_out + band_name
index = band_names[band_name]
x_size = ds.RasterXSize # Raster xsize
y_size = ds.RasterYSize # Raster ysize
driver = gdal.GetDriverByName('GTiff')
arch = driver.Create(filename_output,x_size,y_size,1,gdal.GDT_Float32)
arch.SetGeoTransform(geo_transform)
arch.SetProjection(srs)
arch.GetRasterBand(1).WriteArray(bands_out[:,:,index].astype(np.float32))
del(arch)
print("Band "+band_name+" exported")
In [11]:
rootdir = "Raster/Segmented"
# path to your training data
path_pix = "Entrenamiento/"
# path to your model
path_model = "Data/Models/"
# path to your classification results
path_class = "Data/Class/"
samples = path_pix + "Entrenamiento_raster.tif"
ds_seg1 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B1.tif')
ds_seg2 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B2.tif')
ds_seg3 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B3.tif')
ds_seg4 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B4.tif')
ds_seg5 = gdal.Open('Raster/Segmented/M_2018-05-12_Corte_Seg_B5.tif')
Segmented_image = np.stack([ds_seg1.GetRasterBand(1).ReadAsArray(),
ds_seg2.GetRasterBand(1).ReadAsArray(),
ds_seg3.GetRasterBand(1).ReadAsArray(),
ds_seg4.GetRasterBand(1).ReadAsArray(),
ds_seg5.GetRasterBand(1).ReadAsArray()], axis=2)
Training and classification of the dataset is performed by the next function:
In [12]:
def plot_time(dt):
fig, ax = plt.subplots()
methods = ('RF', 'GBC', 'SVC', 'LSVC', 'KNN')
y_pos = np.arange(len(methods))
performance = dt
ax.barh(y_pos, performance, align='center')
ax.set_yticks(y_pos)
ax.set_yticklabels(methods)
ax.invert_yaxis() # labels read top-to-bottom
ax.set_xlabel('Seconds')
ax.set_title('Performance by method')
plt.show()
# declare a new function
def training(img_ds, samples, path_pix, path_model, path_class):
img = img_ds.copy()
# Load trainig classes
roi_ds = io.imread(samples)
roi = np.array(roi_ds, dtype='int8')
# read in your labels
labels = np.unique(roi[roi > 0])
print('The training data include {n} classes: {classes}'.format(n=labels.size, classes=labels))
# compose your X,Y data (dataset - training data)
X = img[roi > 0, :]
Y = roi[roi > 0]
# Split dataset into trainig an test
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, random_state=42)
# assign class weights (class 1 has the weight 3, etc.)
#weights = {1:3, 2:2, 3:2, 4:2}
weights = None
#weights = "balanced_subsample"
################### Random Forest Classifier ##########################
# build your Random Forest Classifier
# for more information: http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html
ti_rf = time.time()
rf = RandomForestClassifier(class_weight = weights, n_estimators = 500, criterion = 'gini', max_depth = 4,
min_samples_split = 2, min_samples_leaf = 2, max_features = 'auto',
bootstrap = True, oob_score = True, n_jobs = 2, random_state = None, verbose = True)
# now fit your training data with the original dataset
rf = rf.fit(X_train,Y_train)
# Save your Random Forest Model
model = path_model + "model_RF.pkl"
joblib.dump(rf, model)
# call your random forest model
rf = path_model + "model_RF.pkl"
clf = joblib.load(rf)
# Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack)
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :23].reshape(new_shape)
class_prediction = clf.predict(img_as_array)
class_prediction = class_prediction.reshape(img[:, :, 0].shape)
tf_rf = time.time()
dt_rf = tf_rf - ti_rf
#########################Gradient Boost Classifier#############################
# alternatively you may try out a Gradient Boosting Classifier
# It is much less RAM consuming and considers weak training data
ti_gbc = time.time()
GBC = GradientBoostingClassifier(n_estimators = 500, min_samples_leaf = 1, min_samples_split = 4, max_depth = 4,
max_features = 'auto', learning_rate = 0.8, subsample = 1, random_state = None,
warm_start = True)
# now fit your training data with the original dataset
GBC = GBC.fit(X_train,Y_train)
# Save your Gradient Boosting Model
model_GBC = path_model + "model_GBC.pkl"
joblib.dump(GBC, model_GBC)
GBC = path_model + "model_GBC.pkl"
clf_GBC = joblib.load(GBC)
# Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack)
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :23].reshape(new_shape)
class_prediction_GBC = clf_GBC.predict(img_as_array)
class_prediction_GBC = class_prediction_GBC.reshape(img[:, :, 0].shape)
tf_gbc = time.time()
dt_gbc = tf_gbc - ti_gbc
#########################Support Vector Classifier#############################
ti_svc = time.time()
svc = SVC(random_state = None, gamma = 'auto')
# now fit your training data with the original dataset
svc = svc.fit(X_train,Y_train)
# export your SVC model
model_svc = path_model + "model_svc.pkl"
joblib.dump(svc, model_svc)
# call your SVC model
svc = path_model + "model_svc.pkl"
clf_svc = joblib.load(svc)
# Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack)
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :23].reshape(new_shape)
class_prediction_svc = clf_svc.predict(img_as_array)
class_prediction_svc = class_prediction_svc.reshape(img[:, :, 0].shape)
tf_svc = time.time()
dt_svc = tf_svc - ti_svc
#########################Linear Support Vector Classifier#############################
ti_lsvc = time.time()
lsvc = LinearSVC()
# now fit your training data with the original dataset
lsvc = lsvc.fit(X_train,Y_train)
# export your Random Forest / Gradient Boosting Model
model_lsvc = path_model + "model_lsvc.pkl"
joblib.dump(lsvc, model_lsvc)
# call your Linear SVC model
lsvc = path_model + "model_lsvc.pkl"
clf_lsvc = joblib.load(lsvc)
# Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack)
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :23].reshape(new_shape)
class_prediction_lsvc = clf_lsvc.predict(img_as_array)
class_prediction_lsvc = class_prediction_lsvc.reshape(img[:, :, 0].shape)
tf_lsvc = time.time()
dt_lsvc = tf_lsvc - ti_lsvc
######################### KNN Classifier#############################
ti_knn = time.time()
knn = KNeighborsClassifier()
# now fit your training data with the original dataset
knn = knn.fit(X_train,Y_train)
# export your KNN Model
model_knn = path_model + "model_knn.pkl"
joblib.dump(knn, model_knn)
# call your KNN model
knn = path_model + "model_knn.pkl"
clf_knn = joblib.load(knn)
# Classification of array and save as image (23 refers to the number of multitemporal NDVI bands in the stack)
new_shape = (img.shape[0] * img.shape[1], img.shape[2])
img_as_array = img[:, :, :23].reshape(new_shape)
class_prediction_knn = clf_knn.predict(img_as_array)
class_prediction_knn = class_prediction_knn.reshape(img[:, :, 0].shape)
tf_knn = time.time()
dt_knn = tf_knn - ti_knn
########## Plot ROC curves ####################
ax = plt.gca()
rfc_disp = plot_roc_curve(clf,X_test,Y_test, ax=ax)
GBC_disp = plot_roc_curve(clf_GBC,X_test,Y_test, ax=rfc_disp.ax_)
svc_disp = plot_roc_curve(clf_svc,X_test,Y_test, ax=GBC_disp.ax_, alpha=0.8)
lsvc_disp = plot_roc_curve(clf_lsvc,X_test,Y_test, ax=svc_disp.ax_, alpha=0.8)
knn_disp = plot_roc_curve(clf_knn,X_test,Y_test, ax=lsvc_disp.ax_, alpha=0.8)
knn_disp.figure_.suptitle("ROC curve comparison")
plt.show()
ax1 = plt.gca()
ax1.set(xlim=[0,0.7], ylim=[0.4,1])
rfc_disp = plot_roc_curve(clf,X_test,Y_test, ax=ax1)
GBC_disp = plot_roc_curve(clf_GBC,X_test,Y_test, ax=rfc_disp.ax_)
svc_disp = plot_roc_curve(clf_svc,X_test,Y_test, ax=GBC_disp.ax_, alpha=0.8)
lsvc_disp = plot_roc_curve(clf_lsvc,X_test,Y_test, ax=svc_disp.ax_, alpha=0.8)
knn_disp = plot_roc_curve(clf_knn,X_test,Y_test, ax=lsvc_disp.ax_, alpha=0.8)
knn_disp.figure_.suptitle("ROC curve comparison - ZOOM")
plt.show()
plot_time([dt_rf, dt_gbc, dt_svc, dt_lsvc,dt_knn])
return class_prediction, class_prediction_GBC, class_prediction_svc, class_prediction_lsvc, class_prediction_knn
In [13]:
class_prediction, class_prediction_GBC, class_prediction_svc, class_prediction_lsvc, class_prediction_knn = training(Segmented_image, samples, path_pix, path_model, path_class)
The scikit-klearn methods do not support NaN data, so we can't replace pixels with zero values by using np.nan, so when classifying, pixels with a zero value are included in one of the classes, in this case, in the late blight class. To clean the results, we use the mask created previously.
In [14]:
from matplotlib.colors import from_levels_and_colors
cmap, norm = from_levels_and_colors([0,0.5,1,2],['black','red','yellow'])
plt.figure(3, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231) ,plt.imshow(class_prediction, cmap=cmap),plt.title('RFC')
plt.subplot(232) ,plt.imshow(class_prediction_GBC,cmap=cmap),plt.title('GBC')
plt.subplot(233) ,plt.imshow(class_prediction_svc,cmap=cmap),plt.title('SVC')
plt.subplot(234) ,plt.imshow(class_prediction_lsvc,cmap=cmap),plt.title('LSVC')
plt.subplot(235) ,plt.imshow(class_prediction_knn,cmap=cmap),plt.title('KNN')
plt.show()
In [15]:
cmap, norm = from_levels_and_colors([0,0.5,1,2],['white','#EA8A00','green'])
rf_classification = class_prediction*filtered_mask
gbc_classification = class_prediction_GBC*filtered_mask
svc_classification = class_prediction_svc*filtered_mask
lsvc_classification = class_prediction_lsvc*filtered_mask
knn_classification = class_prediction_knn*filtered_mask
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231) ,plt.imshow(rf_classification, cmap=cmap),plt.title('RFC')
plt.subplot(232) ,plt.imshow(gbc_classification,cmap=cmap),plt.title('GBC')
plt.subplot(233) ,plt.imshow(svc_classification,cmap=cmap),plt.title('SVC')
plt.subplot(234) ,plt.imshow(lsvc_classification,cmap=cmap),plt.title('LSVC')
plt.subplot(235) ,plt.imshow(knn_classification,cmap=cmap),plt.title('KNN')
plt.show()
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(231) ,plt.imshow(rf_classification[1400:1700,1200:1400], cmap=cmap),plt.title('RFC Zoom')
plt.subplot(232) ,plt.imshow(gbc_classification[1400:1700,1200:1400],cmap=cmap),plt.title('GBC Zoom')
plt.subplot(233) ,plt.imshow(svc_classification[1400:1700,1200:1400],cmap=cmap),plt.title('SVC Zoom')
plt.subplot(234) ,plt.imshow(lsvc_classification[1400:1700,1200:1400],cmap=cmap),plt.title('LSVC Zoom')
plt.subplot(235) ,plt.imshow(knn_classification[1400:1700,1200:1400],cmap=cmap),plt.title('KNN Zoom')
plt.show()
In [16]:
labels_class = np.unique(class_prediction[class_prediction > 0])
print('The classified data include {n} classes: {classes}'.format(n=labels_class.size, classes=labels_class))
In [17]:
ds_ground_truth = gdal.Open('Validacion/Raster/Ground_truth_2018-05-12_.tif')
# asignacion cada banda
ground_truth_array_raw = ds_ground_truth.GetRasterBand(1).ReadAsArray()
ground_truth_array = ground_truth_array_raw*filtered_mask
srs = ds_ground_truth.GetProjectionRef()
geo_transform = ds_ground_truth.GetGeoTransform()
cmap, norm = from_levels_and_colors([0,0.5,1,2],['white','#EA8A00','green'])
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(111) ,plt.imshow(ground_truth_array, cmap=cmap),plt.title('Ground truth')
plt.show()
plt.figure(1, dpi=300)
plt.subplots_adjust(left=0.0, right=3.0, bottom=0.0, top=3.0)
plt.subplot(111) ,plt.imshow(ground_truth_array[1400:1700,1200:1400], cmap=cmap),plt.title('Ground truth Zoom')
plt.show()
Ground truth details¶
In [18]:
def disease_info(array, geo_transform):
unique, counts = np.unique(array[array!=0], return_counts=True)
pix_width = geo_transform[1]
pix_area = pix_width*pix_width
lb_total_area = counts[0]*pix_area
h_total_area = counts[1]*pix_area
print("Late blight afected area: ", round(lb_total_area,2), "square meters")
print("Healthy potato area: ", round(h_total_area,2), "square meters")
print("Late blight vs Healthy potato area ratio: ", round(lb_total_area/h_total_area,2))
print("Late blight area vs total area: ", round(lb_total_area/(lb_total_area+h_total_area),2))
print("Healthy potato area vs total area: ", round(h_total_area/(lb_total_area+h_total_area),2))
In [19]:
disease_info(ground_truth_array, geo_transform)
Classification accuracy assessment¶
Random forest:¶
In [20]:
y_true = ground_truth_array.ravel()
target_names = ['0', 'Late blight', 'Healthy']
y_pred_rf = rf_classification.ravel()
print(classification_report(y_true, y_pred_rf, target_names=target_names))
Disease details of Random forest¶
In [21]:
disease_info(rf_classification, geo_transform)
Gradient Boost Classifier¶
In [22]:
y_pred_gbc = gbc_classification.ravel()
print(classification_report(y_true, y_pred_gbc, target_names=target_names))
Disease details Gradient Boost Classifier¶
In [23]:
disease_info(gbc_classification, geo_transform)
Support Vector Classifier¶
In [24]:
y_pred_svc = svc_classification.ravel()
print(classification_report(y_true, y_pred_svc, target_names=target_names))
Disease details Support Vector Classifier¶
In [25]:
disease_info(svc_classification, geo_transform)
Linear Support Vector Classifier¶
In [26]:
y_pred_lsvc = lsvc_classification.ravel()
print(classification_report(y_true, y_pred_lsvc, target_names=target_names))
Disease details Linear Support Vector Classifier¶
In [27]:
disease_info(lsvc_classification, geo_transform)
K-Neighbors Classifier¶
In [28]:
y_pred_knn = knn_classification.ravel()
print(classification_report(y_true, y_pred_knn, target_names=target_names))
Disease details K-Neighbors Classifier¶
In [29]:
disease_info(knn_classification, geo_transform)
References¶
Glasbey, C.A., 1993. An Analysis of Histogram-Based Thresholding Algorithms.
Otsu, N., 1979. A Threshold Selection Method from Gray-Level Histograms. IEEE transactions on systems, man, and cybernetics 9, 62–66.arXiv:1011.1669v3.
Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011